Mountains, climate and niche heterogeneity explain global patterns of fern diversity

1.1 Speciation rate estimation

Lineage speciation rates were estimated using a time‐calibrated phylogeny of 3973 fern taxa (Testo & Sundue, 2016), implementing a branch‐specific birth–death model in RevBayes (Höhna et al., 2016). We accounted for incomplete taxon sampling by assigning family‐level sampling proportions, using species richness estimates from the Pteridophyte Phylogeny Group (PPG, 2016). We specified a prior on the number of rate categories at k = 6 and assigned lognormal distributions for speciation and extinction rates with means equal to 4.64 ((number of taxa in phylogeny/root age)/2) and standard deviations set to 0.5. We allowed for rate shifts to occur across the tree and set a uniform distribution on the number of shifts scaled between 0 and 10. This model was sampled with an MCMC on the University of Florida’s HiPerGator 2.0 high‐performance computing cluster for 24,000 generations, with 4000 generations specified as burn‐in. Convergence and mixing were assessed using Tracer 1.7 (Rambaut et al., 2018) to evaluate run convergence and to ensure that all parameters had an effective sample size of > 200. Tip speciation rates were extracted and used for subsequent analyses.

1.2 Spatial analyses

Species occurrence records were downloaded from GBIF with queries filtered to include only georeferenced fern herbarium specimens, resulting in a total of 1,686,374 records (http://GBIF.org“ GBIF.org (22 April 2020) GBIF Occurrence Download https://doi.org/10.15468/dl.amnnxg”.). To improve the quality of these records, we carried out a multi‐step occurrence cleaning process. First, we used the R package ‘CoordinateCleaner’ (Zizka et al., 2019) to flag and remove records that satisfied any of the following criteria: 0,0 coordinates, coordinates in an ocean, coordinates within 5 km of country centroids and capitals, records within 1 km of biodiversity institutions, and records with reversed latitude and longitude values; this removed 809,228 records. After removing these geographical errors, we resolved the taxonomy by passing all remaining records through a taxonomic thesaurus developed from the Checklist of Ferns and Lycophytes of the World (https://worldplants.webarchiv.kit.edu/ferns/), with additional inputs from Tropicos (www.tropicos.org) and the Pteridophyte Collections Consortium data portal (www.pteridoportal.org). The taxonomic thesaurus included a total of 59,448 published names (10,657 accepted names). All records were matched against this list and records with names that were classified as synonyms were updated to the corresponding accepted name; this resulted in name changes for 197,875 records. A total of 2386 records could not be matched to an accepted name and were excluded. Next, we used the R package ‘rbokeh’ (Hafen & Continuum Analytics, 2016) to manually generate interactive maps of species occurrences (visualized at the genus level) to identify further erroneous records, species outside of their native range and flag suspicious records for additional inspection; an additional 29,299 records were manually removed, resulting in a total of 845,878 records representing 7865 total species.

We summarized regional fern diversity by associating cleaned occurrence records with 100 × 100 km grid cells using the ‘join attributes by location’ tool in QGIS 3.10.1 and trimming the resulting data to retain distinct species per grid cell. For each grid cell, we calculated species richness and weighted endemism (Williams & Humphries, 1994), which represents the mean of the inverse range sizes of the taxa represented in each grid cell weighted by species richness. Weighted endemism aims to avoid issues with the use of arbitrary range sizes used to define an endemic taxon (Slatyer et al., 2007). All maps were depicted with an equal‐area (Mollweide) projection.

Hotspots of species richness and endemism were identified using the Gi* local statistic (Getis & Ord, 1996) implemented in QGIS with the Hotspot Analysis plugin (Oxoli et al., 2017). This approach allows for the detection of spatial autocorrelation of a given variable(s) (in this case, species richness and endemism, which were calculated for each grid cell), thereby allowing clusters of user‐defined areas (in this case, 100 × 100 km grid cells) of high values to be identified. We allowed for distance band optimization with a distance step of 10 km and identified global fern diversity hotspots as regions of adjacent grid cells with Z scores ≥ 1.65 for either species richness and endemism or both. Once hotspot cells were identified, we also designated nearby ‘non‐hotspot’ areas in three broad regions for comparative purposes. These non‐hotspots included all grid cells that were not identified as hotspots of either richness or endemism in the following areas: (1) all of continental America from 23ºN to 35ºS (tropical and southern subtropical America) and (2) all of Africa south of 10ºN (sub‐Saharan Africa), and Asia/Malesia from 70–155ºE and 30ºN–10ºS (eastern and southern Asia). Species richness and endemism estimates were obtained by extracting occurrence data within the cells comprising each hotspot; species were considered endemic to a hotspot if no occurrence records for the species were recovered outside of that hotspot. To generate estimates of species richness per elevational band, we extracted elevation values from the GMTED DEM layer (see ‘Section 2.3’) for each occurrence point and then sorted records by that value.

Mean speciation rates of lineage assemblages in each grid cell (average grid cell speciation rates) were calculated based on the 3218 species for which rates were estimated (these were species represented both in our occurrence data and the dated phylogeny of Testo & Sundue, 2016). A total of 604,748 occurrence points for species with speciation rate data were obtained from our total dataset and a single representative per grid cell was retained to calculate average grid cell speciation rates. Using these occurrences and time‐calibrated phylogeny, we also calculated spatial phylogenetic indices of phylogenetic endemism (PE) and phylogenetic diversity (PD) for each grid cell. PD is the sum of the total branch lengths that connect sets of taxa to a most recent common ancestor (Faith, 1992). PE is the range‐size‐weighted equivalent of PD (Laffan et al., 2016; Rosauer et al., 2009). These indices were calculated for use in downstream analyses, which are described below.

We identified centres of neo‐ and palaeo‐endemism with the CANAPE analysis (Mishler et al., 2014), which identifies areas with different endemism signals through a two‐step process of evaluating relative PE (RPE) values. RPE is the ratio comparing the observed PE from the original tree to PE where all branches of the tree are of equal length (PE_alt). In the first step, we used Biodiverse 3.0.0 (Laffan et al., 2010) to calculate PE, PE_alt and RPE values for all grid cells in our dataset and identify grid cells with significantly high values for either PE or PE_alt; significance (two‐tailed test, α = 0.05) was assessed by comparison to a null distribution of RPE values generated with 999 iterations of the ‘rand_structured’ function. In the second step, grid cells with significantly high values of either PE or PE_alt (or both) were then divided into four endemism classes based on their RPE scores. Grid cells with significantly high RPE are considered to be centres of palaeo‐endemism, cells with significantly low RPE are considered to be centres of neo‐endemism, cells with non‐significant RPE but significantly high PE and PE_alt are considered mixed endemism sites; in cases where both PE and PE_alt are significant at the α = 0.001 level, the cell is considered super‐endemic.

1.3 Environmental and climatic data

To compare environmental niche space of diversity hotspots and non‐hotspots, we obtained climate, soil and topographic data from CHELSA (Karger et al., 2017), SoilGrids (Hengl et al., 2017) and the USGS GMTED2010 (Danielson & Gesch, 2011) databases. The variables used were as follows: annual mean temperature, minimum temperature of the coldest month, annual precipitation, precipitation seasonality, soil pH, soil carbon content and mean elevation. Climate and elevation layers used were at a 30‐arc‐second resolution; soil variables were at a 15 cm depth and a 250 m resolution. Values for each variable were extracted at specimen occurrence points using the ‘sample raster values’ function in QGIS; values were then associated with their corresponding 100 × 100 km grid cell using the ‘join attributes by location’ tool. Environmental niche heterogeneity per grid cell was calculated as the mean of the standard deviations of all environmental variables extracted at species occurrences within a given grid cell. We tested for correlations between individual environmental variables and richness and speciation rate using multiple regressions in R (R Core Team, 2018). We initially conducted an ordinary least squares (OLS) regression and then to account for spatial autocorrelation we used a simultaneous autoregressive (SAR) lag model (Kissling & Carl, 2007) implemented using the ‘spdep’ package (Bivand & Piras, 2015). For each region, we used the number of nearest neighbours that optimized the fit of the model based on Akaike’s information criterion (AIC). We compared the SAR versus OLS model with AIC. In the linear regressions, we used the R² (or Nagelkerke pseudo‐R² in the case of the SAR; Nagelkerke, 1991) value from the best fitting model. Finally, we estimated the relative contribution of each environmental variable to global fern species richness using a multiple regression and accounting for spatial autocorrelation.

3.1 Tropical mountains influence global patterns of species richness

Mountains, especially in the tropics, harbour disproportionately high levels of species richness and endemism relative to the land area they occupy (Rahbek, Borregaard, Colwell, et al., 2019). We demonstrate that ferns follow this pattern with striking prominence: 58% of all species occur in eight principally montane tropical and subtropical hotspots which together comprise just 7% of Earth’s land area (Figure 2; Appendix S2). This exceptional richness is mirrored by increased climate space in these regions. For example, the tropical Andes harbour more than 50% of global temperature–precipitation variation, despite covering just 1.2% of Earth’s land surface (Figure 5). Tropical and subtropical mountains with steep elevational gradients lead to greater ecological niche space compared to surrounding lowlands with lower fern diversity, and we recover a scalable, nearly universal, pattern of local fern species richness being positively correlated with available niche space (Figures S1.9 and S1.10). This is compounded by the fact that the volume of mountainous areas is much greater than the planimetric surface area they occupy. The pattern of outsized climate space occupancy is even more pronounced when only montane portions of these areas are considered: hotspot areas above 1000 m in elevation cover < 2.0% of Earth’s land surface but harbour over 45% of Earth’s fern diversity (Figure 2).

Combining all records (i.e., in hotspot and non‐hotspots), we find a strong latitudinal diversity gradient (Figure 3b) and a middle to upper elevational peak in species richness (Figure 3a), as reported in previous transect‐based studies (Hernández‐Rojas et al., 2020; Kessler et al., 2011; Kluge et al., 2006; Weigand et al., 2020). While the mid‐elevation richness peak is a well‐established phenomenon observed in many tropical plant and animal lineages (Rahbek, 1995), the pattern is particularly strong and consistent in ferns (Cardelús et al., 2006; Kessler et al., 2011; Watkins et al., 2006). Contention about the factors underlying the general mid‐elevation richness peak aside (Colwell & Lees, 2000), understanding why ferns adhere to this pattern so strongly is critical to explain variation in the group’s richness at broad scales.

We propose that the exceptional fern species richness in tropical middle to upper elevation forests (Figure 3a) is largely due to the prominence of many epiphytic fern lineages in these regions. Although epiphytic ferns do not appear to diversify more rapidly than terrestrial species (Schuettpelz & Pryer, 2009; Testo & Sundue, 2018), nearly perhumid conditions, due to significant water input from ‘horizontal precipitation’ (Vogelmann, 1973), allow for the persistence of a diverse vascular epiphyte flora in tropical cloud forests (Cardelús et al., 2006; Küper et al., 2004). Taxonomic representation among vascular epiphytes is highly skewed with a few successful lineages such as orchids and bromeliads comprising the majority of diversity (Gentry & Dodson, 1987). Along with orchids, ‘ferns are the predominant vascular epiphytes nearly everywhere’ (Gentry & Dodson, 1987), and can comprise up to 25% of the vascular epiphyte flora in some tropical cloud forests (Cardelús et al., 2006), despite making up less than 3% of global vascular plant diversity (PPG, 2016). In many middle to upper elevation tropical forests, epiphytic ferns are more diverse than their terrestrial counterparts; for instance in the cloud forests of Costa Rica epiphytic species can comprise 84% of local fern diversity (Watkins et al., 2006). The success of ferns in epiphytic niches is also clear when viewed from an evolutionary perspective: several of the most diverse fern clades (e.g. Asplenium, Elaphoglossum, Hymenophyllaceae, Polypodiaceae) are principally epiphytic (Schuettpelz & Pryer, 2009; Testo & Sundue, 2016). Though epiphytic ferns occur across a range of elevations, and even into temperate regions, their exceptional diversity in montane forests is clearly a principal component of the overall prominence of ferns in middle to upper elevation sites in the tropics.

3.2 Geographical patterns and drivers of speciation

On a global scale, we find that rapidly speciating lineages are primarily found in middle to lower latitudes (Figure S1.6), particularly in montane sites in East Asia and northern South America (Figure S1.5). Speciation rates of lineage assemblages in most regions are higher within hotspot grid cells compared to those in adjacent non‐hotspot grid cells (Figure S1.8). Two hotspots harbouring some of the most rapidly speciating lineages are the tropical Andes and East Asia (Figures S1.5 and S1.7). These sites correspond closely to mountain systems that have undergone recent and rapid orogeny and include some of the steepest elevational and climatic gradients on Earth (Clark et al., 2005; Hoorn et al., 2010). In addition to ferns (Figure S1.5), these same mountains also host other rapidly diversifying lineages (Madriñán et al., 2013; Xing & Ree, 2017). Recent studies have demonstrated a correlation between rates of mountain uplift and the exceptional speciation in plant groups such as bellflowers (Lagomarsino et al., 2016), bamboos (Ye et al., 2019), clubmosses (Testo et al., 2019) and orchids (Pérez‐Escobar et al., 2017). We find the tropical Andes and Hengduan Mountains include some of the youngest and most rapidly speciating fern lineages, such as Serpocaulon and Elaphoglossum in the Andes, and Athyrium and Lepisorus in the Hengduan Mountains (Appendix S2). However, we also find examples of rapid speciation outside of our circumscribed hotspots: both the Cape Region of Southern Africa and parts of Australia harbour rapidly diversifying fern assemblages (Figure S1.5). In these cases, it appears that, although there is evidence of recent radiations, overall fern diversity remains modest (Figure 1).

Within hotspots, we found a strong relationship between speciation rates of lineage assemblages and environmental niche heterogeneity (Figure 6); this is most evident in the tropical Andes, Mesoamerica and East Asia. Globally, we find sites of rapid speciation are generally in topographically complex regions (Figure S1.5). Geological, geographical and climatic factors underlie processes that play an important role in promoting speciation, including orogeny‐induced generation of novel niche space, rapid climatic turnover along steep elevational gradients, and local variation in soil and bedrock types (Condamine et al., 2018; Givnish et al., 2014; Hoorn et al., 2010, 2013, 2018; Rahbek, Borregaard, Antonelli, et al., 2019).

3.3 Factors influencing patterns of endemism

Both endemism and species richness are higher in montane sites when considered on a per‐unit‐area basis (Figure 2). However, total numbers of endemic species in most hotspots were similar in montane and lowland habitats (Figure 2); this contrasts sharply with a general pattern of increased endemism at high elevations that has been reported for many groups of organisms (Rahbek, Borregaard, Colwell, et al., 2019). Reasons for this departure remain uncertain, but we propose several possibilities. First, the high dispersibility of fern spores may allow species inhabiting patchily distributed montane habitats to maintain sufficient gene flow between populations to avoid isolation on ‘sky islands’ (Barrington, 1993), which is thought to be a principal mechanism by which narrowly endemic species form in other montane groups (Hughes & Eastwood, 2006). Second, available evidence suggests that fern lineages tend to show strong affinity to soil types, and species distributions at local scales may largely be mediated by edaphic features (Lehtonen et al., 2017; Tryon, 1972). Because diversity of soil and bedrock types is largely independent of elevation at broad spatial scales (Antonelli et al., 2018; Rahbek, Borregaard, Antonelli, et al., 2019), edaphic constraints on species distributions may help explain the patterns observed here. Historical climate factors not evaluated in our analyses may have played an important role in shaping patterns of endemism along elevational gradients. As hypothesized by Kluge et al. (2006) for Mesoamerican ferns, upwards shifts of forest assemblages during periods of global warming in the Quaternary may have reduced the number of high‐elevation endemics through local extinctions on isolated mountaintops.

Although levels of endemism are high across global hotspots (Figure S1.2), the type of PE varies by site (Figure S1.4). Parts of the tropical Andes, Madagascar and East Asia include sites with high levels of neo‐endemism (Figure S1.4); Malesia, in contrast, harbours most of the world’s centres of palaeo‐endemism and mixed‐endemism (areas with both high neo‐endemism and palaeo‐endemism; Figure S1.4). The distribution of these endemism centres provides insight into the disparate evolutionary histories of ferns in these regions. Within the tropical Andes and the Hengduan Mountains, evidence suggests that areas with incidence of neo‐endemism are due to recent in situ speciation at high elevations (Figures S1.5 and S1.7; Sánchez‐Baracaldo & Thomas, 2014; Wei et al., 2018; Xing & Ree, 2017). As a rule, these habitats are geologically complex and distribution of their vegetation is strongly impacted by climatic variation along an elevation cline (Clark et al., 2005; Hoorn et al., 2010). Dynamics between climatic variation across, and seasonal stability within, elevational transects in these tropical and subtropical mountains is hypothesized to reduce the range size of many groups of organisms leading to high allopatric speciation (Janzen, 1967; Polato et al., 2018). Further study is needed to better understand how ‘sky islands’ and climate‐driven shifts in their extent impacts patterns of speciation and endemism in ferns and other groups; methods that incorporate spatially explicit models of habitat connectivity through time (Flantua et al., 2019) may prove to be particularly insightful.

In contrast, the predominance of palaeo‐endemism in Malesia highlights the importance of isolation and long‐term climatic stability in maintaining the diverse and unique species assemblages found in this hotspot (Figure S1.4). Although we find evidence of rapid, recent speciation in some areas (e.g. the Papuan Highlands and Borneo; Figure S1.5), Malesia has the lowest mean speciation rate of any of our hotspots (Figure S1.7). Compared to other global hotspots, the high fern diversity observed in Malesia appears to be largely due to the persistence of ancient endemic lineages. Most Malesian ferns are found nowhere else on Earth (Figure 2), including members of relictual, once widespread, families such as Cystodiaceae and Matoniaceae, which are represented in the fossil record dating back hundreds of millions of years (Choo & Escapa, 2018; Regalado et al., 2018). The importance of Malesia as both a biodiversity hotspot and ‘museum’ has been noted for other groups (Steenis, 1979) and reflects the signature of the region’s complex palaeogeographical and climatic histories on its biota. Since the Cretaceous, Malesia has experienced a complex history of land submersion and emergence events, resulting in a continually changing mosaic island chains with varying levels of connectivity (Lomolino et al., 2006). This highly dynamic landscape history led to a multitude of local speciation, extinction and migration events that resulted in conspicuous patterns of distributions observed today, and have contributed to exceptional diversification bursts in numerous lineages (Moyle et al., 2016; Thomas et al., 2012). At the same time, Malesia has acted as a refuge for many relictual lineages that were formerly widespread, suggesting that this region has experienced reduced rates of extinction compared to continental hotspots. The persistence of ancient groups on these and other islands (such as Madagascar) is likely aided by the thermal stability of insular systems and lower rates of biotic turnover (Cronk, 1997; Fjeldså et al., 2012).

3.4 Conclusions, caveats and future directions

Recent work focusing on tropical mountains has highlighted the geologic and climatic factors that help generate and maintain their exceptional biodiversity, but most broad‐scale analyses have focused on a few well‐known groups, such as mammals, birds and amphibians (Antonelli et al., 2018; Perrigo et al., 2020; Quintero & Jetz, 2018; Rahbek, Borregaard, Antonelli, et al., 2019). In this study, we close a significant gap in our knowledge with a global analysis of patterns of diversity in a major plant lineage: ferns. By integrating phylogenetic, ecological and georeferenced data, we demonstrate the following points. By our estimates, 58% of global fern diversity is concentrated in only eight biodiversity hotspots and nearly half (~45%) is concentrated in montane sites (Figures 1 and 2). Second, these hotspots occupy a disproportionate amount of climatic and environmental niche space compared to their total land area, which acts as a strong predictor of its species richness (Figure 5; Figures S1.9 and S1.10). Third, species assemblages with the highest rates of speciation are generally found at middle to lower latitudes and higher elevations (Figure S1.6), especially in areas with high environmental heterogeneity (Figure 6). Finally, patterns of endemism vary across hotspots: neo‐endemism is prevalent in areas with high, young mountains, whereas tropical islands with stable climates and relatively low biotic turnover, act as refugia for ancient lineages (Figure 3). Taken together, these findings suggest that mountain building and associated formation of heterogeneous environments promote rapid species proliferation through reduced dispersal and gene flow, while long‐term ecological and climatic stability helps maintain global diversity in this ancient clade by buffering against extinction in formerly widespread groups.

There are several limitations to the present study. First, it is well known that quality of available species occurrence data is highly variable and that misidentified records or those with problematic localities can bias biodiversity analyses (Zizka et al., 2019). Although we made extensive efforts to ensure the quality of the species occurrence data used through a series of automated and manual cleaning steps—resulting in the removal of 50% of the records in the original dataset—some problematic records almost certainly remain in our final dataset. In addition, political and geographical biases in data availability are evident. For example, few records of ferns from India are available on GBIF, despite the fact that its fern flora is diverse and relatively well‐studied (Fraser‐Jenkins et al., 2017). Similarly, the marked difference in record density and estimated species richness between the Indonesian and Papuan sides of New Guinea appear to reflect differences in collection effort and digitization, rather than true diversity patterns (Cámara‐Leret et al., 2020). Unsurprisingly, collection effort and inferred species richness are positively correlated globally (Figure S1.1); nonetheless, our estimates of species richness are similar to species counts obtained from floristic treatments for most regions (Appendix S2). These challenges are shared by other studies relying on global species occurrence datasets (Yesson et al., 2007) and interpretation of our findings should be made with appropriate caution. Our findings are also shaped in part by the phylogeny used in several analyses, ranging from speciation rate estimation to characterization of phylogenetic endemism centres. While this phylogeny enabled us to maximize the number of taxa with evolutionary data, more than half of known fern species are yet to be placed in a phylogeny, and some genera remain unsampled altogether (Appendix S2). In addition, although the phylogeny used in this study is largely topologically concordant with those of other studies (Lehtonen et al., 2017; PPG, 2016; Schuettpelz & Pryer, 2009), uncertainty in relationships among lineages and with respect to divergence times could impact some of the findings reported here. Together with improved species occurrence data, phylogenies with improved taxon sampling, increased numbers of loci and more sophisticated divergence time estimation methodologies should provide improved resolution of the patterns explored here.

With these caveats considered, we suspect that our main conclusions are valid. The hotspots of richness and endemism that we recover correspond closely to those reported previously in floristic studies as well as plot‐based analyses examining fern diversity at local to global scales (Hernández‐Rojas et al., 2020; Kessler et al., 2011; Kluge et al., 2006; Moran, 1995, 2008; Nagalingum et al., 2015; Tryon, 1972; Weigand et al., 2020; Xing & Ree, 2017). In addition, our findings are similar to those found for other taxonomic groups and with varying methodological approaches (Antonelli et al., 2018; Myers et al., 2000; Perrigo et al., 2020; Rahbek, Borregaard, Antonelli, et al., 2019). In light of this, while we anticipate that insight into global fern diversity will improve as more data become available for currently underrepresented regions, we see little reason to expect significant departures from the broad‐scale patterns reported here.

We anticipate that future work will continue to investigate drivers of species diversity by leveraging methodological advances and more robust datasets. In particular, generating a comprehensive understanding of the importance of mountains for biodiversity requires improved integration of biology and geology. Though our findings add to a substantial body of work highlighting the exceptional richness, endemism and rapid diversification of montane groups of organisms, we did not attempt to assess how the timing and rates of mountain uplift events correspond with patterns of lineage diversification. Environment‐dependent models of diversification allow for such correlations to be tested (Condamine et al., 2018), but robust estimates of uplift rates are lacking for many mountain systems and applying these models to a global analysis remains challenging. Similar models can be extended to better integrate species diversification estimates with other important historical factors, including changes in global climate and shifts in the presence and connectivity of islands, both oceanic islands and high‐elevation ‘sky islands’. Integrating these models with improved species‐level data will prove particularly valuable as we seek to better understand what factors influence species distributions, patterns of niche evolution and lineages’ evolutionary trajectories through time.